Refinement and Coarsening of Parameterization for the Estimation of Hydraulic Transmissivity

This paper deals with parameterization of hydraulic transmissivity during its estimation. Estimating transmissivity is based on the minimization of an misfit function defined as a leastsquares misfit between measured data and the model output. The transmissivity is assumed to be a piecewise constant space dependent function. Refinement indicators, indicating the effect on the optimal data misfit of adding degrees of freedom to a current set of parameters, are calculated. Theses indicators lead to know where the hydraulic transmissivity needs to be made discontinuous and thus allows for the introduction of degrees of freedom one at a time, preventing from overparameterization. INTRODUCTION When estimating hydraulic transmissivities, which are space dependent coefficients in a parabolic partial differential equations, we minimize a misfit function defined as a least-squares misfit between measurements and the corresponding quantities computed with a chosen parameterization of the transmissivities. One of the difficulties in solving this problem is that, because of the high cost of experimental measurements, the data is usually insufficient to estimate the value of the hydraulic Address all correspondence to this author. transmissivity in each cell of the computational grid. Therefore we have to find a parameterization of the transmissivity which reduces the number of unknowns. We refer to (Sun, 1994; Eppstein-Dougherty, 1996) for a presentation of the parameterizations which are the most commonly used in hydrogeology. Lately multiscale parameterizations (Liu, 1993; ChaventLiu, 1989) have provided a first answer to the problem of choosing the discretization of distributed parameters. With such an approach the parameter estimation problem is solved through successive approximations by refining the scale and the process is stopped when the refinement of the scale does not induce a significant decrease of the misfit function. This method has already brought interesting results in various problems of parameter estimation (Chardaire-Rivière et al., 1990; Chardighy et al., 1996). However, when going from the current scale to a finer one, degrees of freedom are added uniformly in the domain of calculation, and so this approach can lead to overparameterization in the case where only local refinements are needed. The method using refinement and coarsening indicators, that we present in this paper, will avoid such a drawback. At a given refinement level the parameters are estimated by minimizing the least-squares misfit to the data. Then we compute refinement and coarsening indicators which indicate the effect on the optimal data misfit of adding or removing some degrees of freedom. A 1 Copyright  1999 by ASME variant of this method has already been presented in (ChaventBissel, 1998). We apply this technique to the estimation of the distributed transmissivity parameter T in the partial differential equation S ∂Φ ∂t +div ( T grad Φ) = Q dans Ω; (1) where Ω is a domain in R2 . This equation governs a twodimensional groundwater flow in an isotropic and confined aquifer, subject to initial and boundary conditions Φ = Φd on ΓD; ( T grad Φ) n = qd on ΓN ; Φ(0) = Φ0 in Ω; (2) where ΓD and ΓN is a partition of the boundary of Ω supporting respectively Dirichlet and Neumann conditions and where Φ = piezometric head, S(x;y) = storage coefficient, T (x;y) = hydraulic transmissivity; Q(x;y; t) = distributed source term, Φd ;qd given boundary piezometric head and source terms, Φ0 = given nitial piezometric head, n is the unit normal vector to Ω: The domain Ω is discretized with a mesh T and equations (1),(2) are approximated by a mixed-hybrid finite element method, which is suitable to the case where the parameter T have discontinuities (Chavent-Roberts, 1997). Our problem is to search for a piecewise constant hydraulic transmissivity T . Both the zonation (the partition of Ω into zones where T is constant) and the value of T on each zone are to be determined and our aim is to be able to explain the data with a number of zones as small as possible. We define our misfit function by: J(T ) = 1 2 ∑i; j jΦ(x j; ti) Φi; jj 2 : (3) where Φi; j is the piezometric head measured in the point x j at the time ti and Φ(x j ; ti) is the model output for current transmissivity values. In order to have a computational cost independant of the number of transmissivity values to estimate, we use a Gaussnewton algorithm which is known to be an efficient optimization method in the case of large number of unknown parameters (Lemarechal et al., 1997). The gradient of J is computed by the adjoint state method (Sun, 1994). Actual implementation of the calculation of the the piezometric head, the misfit function and its gradient has been achieved through automatic program generation (Jegou, 1997). 1 REFINEMENT AND COARSENING INDICATORS Usually the data are in insufficient quantity to estimate the unknown parameter in each cell of the computing mesh T . To every partition (zonation) Z of Ω, we associate the space of hydraulic transmissivities T which are constant on each zone of Z. In practice we suppose that the boundary of a zone in Z is made of edges or diagonals of cells of T . The degrees of freedom of the transmissivity to be estimated are the values of the transmissivity in each zone of Z. We want this number to be significantly lower than the number of elements in T . Our goal is that the refinement and coarsening indicators technique will allow us to construct Z with the lowest possible number of zones. Given a current parameterization partition, these indicators will tell us where to insert new discontinuities of the transmissivity and which ones can be removed. 1.1 Refinement indicators We shall describe refinement indicators on an example. Let (P1) be an initial problem where the hydraulic transmissivity is constant in all the domain (figure 1(b)). So we have only one value of the transmissivity to estimate, which is done by minimizing the one variable misfit function J. We note (P2) the parameter estimation problem where we consider that the domain is made of two zones Z1 and Z2 having different transmissivities (figure 1(c)). We note T = (T 1 ;T 2 ) the solution of (P2) obtained by minimizing the corresponding two variable misfit function. If B=T 1 T 2 were known, then the solution of minimizing J under the constraints T1 T2 = B is necessarily T = (T 1 ;T 2 ), the solution of (P2) and, if B = 0 then T = T 1 = T 2 is the solution of (P1). In general if we introduce n new transmissivity values to es(c) (a) (b) Computing mesh T Partition Z1 associated to problem (P1) Z Partition Z2 associated to problem (P2) Z1 Z2 Figure 1. An example of parameterization refinement timate, the constraint condition T1 T2 = B is generalized to 2 Copyright  1999 by ASME AT = B where B is the discontinuity vector (σi;i+1)1 i n and A is an n (n+1)matrix defined by A = 2 6664 1 1 0 0 . . . 0 1 1 . . . . 0 0 1 1 . . . . . . . . . . 0 0 . . . 1 1 3 7775 : To the problem of minimizing the misfit function under these constraints J(T ) = min AT=B J(T ), we associate the Lagrangian function defined by LB(T;λ) = J(T )+< λ;AT B > (4) where λ is the Lagrange multiplier associated to the constraint AT =B. Then the Lagrange condition ensures that T is obtained by solving ∂LB ∂T (T ;λ ) = ∇J(T )+A λ = 0; ∂LB ∂λ (T ;λ ) = AT B= 0: (5) If we denote by J B = J(T ) =LB(T ;λ ) the optimal misfit associated to the right hand side B of the constraint, we deduce from (4) and (5) that ∂J ∂B jB=0 λ = ∂L ∂B (T ;λ )jB=0 = 0: (6) Therefore the Lagrange multiplier gives us the sensitivity of the optimal data misfit J B to the perturbation B. For this reason we call λ a refinement indicator. It can be easily deduced from equation (5). In the case of our example, to the refinement shown in figure 1(c) we associate an indicator λ = ∂J ∂T1 (T ) = ∂J ∂T2 (T ): (7) And without solving (P2), this indicator indicates if the suggested refinement in figure 1(c) is likely to induce an important decrease of the misfit function. A particular case is when the suggested refinement consists in dividing the domain in four zones of different transmissivities like in fig. 2. Discontinuities between values of transmissivity in different zones can be writen: σ1;2 = T1 T2; σ2;4 = T2 T4; σ3;4 = T3 T4; σ1;3 = T1 T3: Z1 Z2